Cooling curves for neutron stars with hadronic matter and quark matter 
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The thermal evolution of isothermal neutron stars is studied with matter both in the hadronic 
phase as well as in the mixed phase of hadronic matter and strange quark matter. In our models, 
the dominant early-stage cooling process is neutrino emission via the direct Urea process. As a 
consequence, the cooling curves fall too fast compared to observations. However, when superfluidity 
is included, the cooling of the neutron stars is significantly slowed down. Furthermore, we find that 
the cooling curves are not very sensitive to the precise details of the mixing between the hadronic 
phase and the quark phase and also of the pairing that leads to superfluidity. 
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I. INTRODUCTION 

Neutron stars arc natural laboratories for physics un- 
der extreme conditions with their extremely high densi- 
ties, powerful energy emission, large magnetic fields, and 
millisecond rotation periods. At the densities near the 
surface of such a star, atoms break apart into nuclei and 
electrons. At higher densities, the electrons neutralize 
with the protons in the nuclei to form neutrons. These 
stars thus consist of a large fraction of neutrons and arc 
supported from gravitational collapse by the neutron de- 
generacy pressure, from which the neutron star derives 
its name. 

However, at high densities the existence of more exotic 
particles is expected. These particles are generated by 
processes which produce strangeness, such as 



A" + 



(1) 



where n is the neutron. A" the Lambda hyperon, and 
the strange meson. The strange meson can decay 
via various weak processes and the final products, usu- 
ally photons and neutrinos, leak out of the star. There- 
fore, the reverse process is reduced and some net amount 
of strangeness survives in the dense core of the star [l| . 
It is generally believed that these strangeness carrying 
particles, called hyperons, can exist in the center of neu- 
tron stars. They coexist with the nucleons as well as the 
leptons e~ and ^~ . The interactions between them are 
dominated by the complicated nuclear force whose car- 
riers are the mesons. We refer to such a system as the 
hadronic phase of matter. The system with only nucle- 
ons and leptons, i.e., without hyperons, is referred to as 
the nuclear phase of matter. 

At even higher densities, as a consequence of asymp- 
totic freedom, quarks become deconfined from the 
hadrons. Therefore, strange quark matter, which con- 
sists oi u, d and s quarks, may also exist in the neutron 
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FIG. 1: The schematic QCD phase diagram as a function of 
temperature T and baryon chemical potential /is. The solid 
lines denote flrst- and second-order phase transitions, whereas 
the dashed lines denote smooth crossovers. According to the 
typical temperature of a proto-neutron star of T ~ 10^*^ K 
~ 1 MeV, the state of matter in the neutron star should be 
close to the bottom line, as indicated by the grey thick line. 



star core. There may also be a phase-separated mix- 
ture of strange quark matter and hadronic matter in a 
certain range of densities 0, which we call the mixed 
phase. At present, there is still a lot unknown about 
the deconfinement phase transition of quarks. The con- 
temporary knowledge on matter at high densities and 
temperatures is shown by the schematic QCD phase di- 
agram in Fig. [T] Since the baryon chemical potential /is 
is a monotonously increasing function from the surface 
to the center of the neutron star, the ^lb axis can be 
mapped to the stellar radius and the different phases in 
the neutron star are explicitly indicated. 

Because of the limited knowledge on the state of mat- 
ter at high densities and the complexity of the interac- 
tions between the particles, many effective models for 
matter inside neutron stars have been constructed. In 
general, these models supply an equation of state, which 
determines the particle composition of the neutron star. 
The observation of neutron stars, in turn, constrain these 
models. For example, if an equation of state is too soft 
it is incapable of supporting a very large stellar mass. 
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such that some models become disfavored whenever the 
data on the heaviest neutron star is updated. The stel- 
lar mass is not the only constraint on the models. The 
cooling of neutron stars can be studied from their lu- 
minosity as a function of time. A proto-neutron star is 
born with a typical temperature larger than 10^'' K, after 
which it mainly loses its energy by two processes, namely 
by neutrino emission everywhere inside the star and by 
photon radiation at the surface. In the early stages of 
the thermal evolution of the star, neutrino emission is 
the dominant cooling effect after which photon radiation 
ultimately takes over Q . Since neutrino emission occurs 
everywhere in the neutron star, it provides a probe for 
studying the state of matter inside the star. 

The most efficient neutrino-emission process is called 
the direct Urea (DUrca) process 

n-^p + e^ + De, p + e^^n + i/e- (2) 

This process is only possible if the proton fraction is more 
than a certain threshold in order for it to satisfy energy 
and momentum conservation Q . Historically, the proton 
abundance in neutron stars was underestimated. In that 
scenario, it is reasonable to take into account also the 
modified Urea (MUrca) process, where a bystander helps 
the momentum conservation, for example 

n + n ^ p + n + + i>e, p + n + ^ 71 + n + Ve, 
n+p^p + p + e^+De, P + P + e^ — > ?i + p + t^e- (3) 

Some other processes can have a neutrino emissivity 
which is smaller or comparable to the modified Urea pro- 
cess, such as neutrino bremsstrahlung and plasmon decay 
0. However, all of them are negligible whenever the di- 
rect Urea channel is open. In the hadronic phase, the 
direct Urea processes are quite rich and can be summa- 
rized as 

bi ^ b2 + + i^i, b2 + r^bi+ui, (4) 

where bi and 62 denote two different baryons and rep- 
resents one of the leptons, e~ or fi~ . 

For strange quark matter, the direct Urea processes 
are 

d — > u + + z?i, u + ^ d + vi, 

s^u + l^+Di, u + l^^s + ui, (5) 

which are simply the direct Urea processes of the baryons 
at the quark level. Inside the neutron star, there is no 
threshold for the direct Urea process in quark matter, 
since the Fermi momenta of the two quarks and lepton 
can always satisfy momentum conservation. As men- 
tioned, the neutrino emissivity depends strongly on the 
type of matter contributing to the process, such that the 
thermal evolution of a star directly probes its composi- 
tion. The main goal of this paper is to study the cooling 
of neutron stars with different types of equations of state 
and try to constrain them by comparing with observa- 
tional data. 



A large amount of research on the cooling process of 
neutron stars has been carried out in the last several 
decades, with a strong focus on the nuclear phase of 
matter. To the best of our knowledge, a unifying model 
covering the low-density nuclear phase to the deconfined 
quark phase, consistent with the QCD phase diagram, 
has been less thoroughly explored. The strange quark 
star composed of strange matter, either with or without 
the nuclear crust, has been considered in some papers 
0, Q. However, the cooling behavior of neutron stars 
containing also all hyperons and possibly also a mixed 
phase of strange quark matter and hadronic matter, has 
not been extensively studied. Two possible reasons for 
this can be given. First, the hadronic and mixed equa- 
tions of state are rather soft and are continuously being 
challenged by new data on heavy neutron stars, such as 
PSR J1903-H0327 which has M = 1.67±0.01 Mq and 
PSR J1614-2230 which even has M = 1.97 ± 0.04 Mg 
[3 . Second, the existence of a mixed phase is being ques- 
tioned in view of screening and surface effects [1]. With 
respect to the first concern, we note that most of the ob- 
served neutron star masses still lie below the maximum 
mass a hadronic or mixed equation of state can allow for. 
According to Ref. [lo| , most nearby young neutron stars 
[m have masses no bigger than 1.4 Mq [12|. Therefore, 
at least for the study of the thermal evolution of these 
stars, the hadronic or mixed phase can still be of great 
importance. In fact, because of the uncertainties in the 
interaction between particles at extremely high density, 
it is hard to exactly determine the equation of state at 
high density. Although the equations of state used in this 
paper cannot support neutron stars as massive as those 
reported above, we can still use them for a discussion on 
medium-mass neutron stars with more complicated com- 
positions. As for the stability of the mixed phase, the 
arguments are still indecisive. For example, many de- 
tails of the surface tension, which strongly influences the 
stability calculation, arc still uncertain. Although it is 
expected that screening and surface effects diminish the 
mixed-phase regime, it is far from certain that its exis- 
tence can be excluded [l^ . The mixed phase can in par- 
ticular have significant effects on the cooling behavior, 
especially for the heat transport inside the star. How- 
ever, we show below that its effect will be less important 
after the star has become isothermal, when the thermal 
evolution is determined by the heat capacity and neu- 
trino emissivity integrated over the whole volume, and 
the inner thermal conductivity no longer plays a role. 

The temperature of a neutron star is generally much 
smaller than the typical Fermi energy as a consequence 
of the very high densities in the star. Therefore, super- 
fiuidity may play an important role. According to BCS 
theory, fermions can form Cooper pairs at low tempera- 
tures via an attractive interaction and thereby lower the 
energy of the system. The resulting pairs, which obey 
Bose statistics, can form a Bose-Einstein condensate and 
the system becomes superfluid. Pair formation changes 
the single-particle dispersion around the Fermi surface 
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and consequently the heat capacity and neutrino emis- 
sivity will be influenced. We find that, without superflu- 
idity, the cooling of neutron stars is too fast compared 
with observations. However, by including the effects of 
superfluidity we obtain a more realistic cooling behavior. 

The quark phase is usually referred to as an exotic 
phase of extremely dense matter, in contrast to the nu- 
clear phase or the hadronic phase. Other exotic phases 
have also been proposed, among which the pion and kaon 
condensates have attracted much attention (l4l . [T5j . How- 
ever, the existence of such phases inside the neutron star 
is still an open question. On the one hand, with such con- 
densates, the equation of state is further softened and 
thus the corresponding maximum star mass is reduced 
p^ . which makes such phases less favored when com- 
pared to the observational data of massive stars, as men- 
tioned above. On the other hand, for the cooling process, 
it was reported that a meson condensate can increase the 
neutrino emissivity over the typical modified Urea emis- 
sivity by several orders of magnitude but it is still 
much less efficient than the direct Urea process. Since in 
our calculation the direct Urea process is always present, 
such enhancement from the meson condensate has a neg- 
ligible effect. Besides, the kaon condensation may even 
reduce the pressure and cause the star to collapse into 
a black hole [l^. Therefore, considering all the above 
arguments, we do not include such meson condensates in 
this paper. 

The paper is organized as follows, we first introduce 
in Sec. H our theoretical framework which includes the 
relativistic stellar structure, the unified mean-field model 
describing the state of matter inside the star, and the 
thermal evolution equations. In Sec. HI, we then present 
the numerical results for the cooling curves of neutron 
stars with the nuclear, hadronic, and mixed equations of 
state. Due to the efficiency of the direct Urea process the 
cooling is seen to be too fast compared to observations. 
Therefore superfluidity is included in Sec. IV and, as a 
consequence, the cooling is slowed down and we find a 
much better agreement with observations. 



II. THEORETICAL FRAMEWORK 

At the high densities of importance to neutron stars 
the neutrons, protons, and electrons can be considered 
to be highly degenerate. In particular, the typical Fermi 
temperature of the nucleons is about 10"'^^ K while the 
temperature of a proto-neutron star is only 10^" K. This 
provides a great theoretical advantage since it makes it 
possible to decouple the calculation of the stellar struc- 
ture, the nuclear model, and the thermal evolution of the 
star. In other words, the thermal evolution of the star 
depends on the stellar structure, which in turn is con- 
structed from the nuclear model. The nuclear model will 
be solved in the zero-temperature limit to give the equa- 
tion of state. The stellar structure is then obtained from 
Einstein's equations using this equation of state. 



Although many neutron stars have fast rotation rates 
and strong magnetic fields, due to the complexity and 
richness of these phenomenon they are beyond our 
present discussion. In fact, the effect of rotation should 
be quite small for most of the neutron stars except for 
millisecond pulsars. Therefore, Einstein's equations are 
solved for static stars without magnetic fields. Further- 
more, the star interior is approximated by a perfect fiuid, 
i.e., it is hydrostatic, which is valid as long as the energy 
transported by heat conduction is negligible compared 
to the total energy. Subsequently, to discuss the thermal 
evolution, some properties of the stellar matter, such as 
the heat capacities, conductivities, and emissivities, can 
be obtained by considering a perturbation of the Fermi 
surfaces of the various particles. In the following, we 
deal with these aspects of our theoretical framework sep- 
arately. 



A. The general relativistic profile of a compact star 

The space-time metric for a non-rotating, spherically 
symmetric star can be written as 



e2*('-)c2dt2 - e'^'^^^dr^ - r^idO^ + sin^ Od^b''), (6) 



where $(r) is the metric function related to the gravi- 
tational redshift, exp[— 2A(r)] = 1 — 2Gra{r)/(?r with 
m{r) = 4:TTr''^ p{r')dr' / the gravitational mass en- 
closed within the sphere of radius r, G is the gravitational 
constant, and c is the speed of light. As a consequence 
of the above metric, Einstein's equations for the case of 
a non-rotating, hydrostatic, and spherically symmetric 
star reduce to the Tolman-Oppenheimer-Volkoff (TOV) 
equations, namely 



d$(r) 

dr 
dp{r) 

dr 



1 



dp{r) 



p{r) 



p{r) dr 
Vp{r)]G[m{r) 



- 47rr'^p(r) /c^ 



1 



2Gm(r) 
2 



(7) 



(8) 



where p(r) and p{r) are the pressure and energy density 
at radius r, respectively. Given an equation of state 
the TOV equations can be numerically integrated to pro- 
vide the structure of the star. The integration starts at 
r = with a given pressure p(0) = Pc until p{R) = is 
reached, which defines the radius R of the star. Outside 
the star, r > R, the metric reduces to the Schwarzschild 
form, exp [2$(r)] = exp [-2A(r)] = 1 - 2GM/c^r, where 
AI = m{R) is the total gravitational mass of the star. 
Note that with the help of a local Lorenz transformation 
we are always able to apply a locally inertial coordinate 
system [l3 such that the calculation of the equation of 
state is performed in a homogeneously flat background. 
Nevertheless, the TOV equations take into account the 
effects of general relativity on the structure of the star. 
However, to calculate this structure we need to supply a 
realistic equation of state for the matter inside a neutron 
star. 
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B. Effective model for matter in neutron stars 



rise to the following field equations: 



To obtain an equation of state for stellar matter, we 
apply an effective nuclear model which includes all rel- 
evant kinds of baryons. The interactions between the 
baryons are mediated by three types of effective meson 
fields, namely the scalar meson cr, the vector meson 
and the isovector meson p^. The scalar field is coupled 
to the derivatives of the baryon field as shown in the 
following Lagrangian density : 



E 



he 



1 + 



b 



= 0, 



he 



*fc-0, (10) 

(11) 

(12) 
(13) 



^6 ihc^^d^ - guoblp^^^ 



hdfj_ad^a 
2? 



fiw^^w^'' mlujf.uj^' hpf.^p^" 



+ 



2hc 4c3 2hc 4c3 

Y.Mificj^d^' -mic')^lji, (9) 
I 



2hc 



where cOq = d/dt, the summation over b and I is over 
all contributing baryon fields ipb and lepton fields Vii 



and in the meson kinetic terms = df^uji, 



dyUJfi and 



p^n = d^pn — di^Pfj, are the antisymmetric field strengths. 
The derivative coupling of the a field was first intro- 
duced by Zimanyi and Moszkowski [isj to remove the 
problem of a too small or even negative reduced baryon 
mass (commonly referred to as the effective mass, how- 
ever, this name is used in this paper only for the effective 
fermion mass on the Fermi surface) at high density in the 
standard Walecka model [l^ . 

This model can be solved within the mean-field ap- 
proximation. In this approximation all the meson fields 
are replaced by their ground-state expectation values, 
which are constants in space-time and their spatial com- 
ponents are zero because the system is assumed to be 
homogeneous and isotropic. Also the charged compo- 
nents of the isospin vector p^ and p~ , whose sources are 
the off-diagonal currents of the baryon fields, are zero. 
Therefore, only the three constant fields cr, w° and p§ 
survive in the Euler-Lagrange equations. In the follow- 
ing, we simply denote the last two as w and p. Further- 
more, the constant scalar field a can be absorbed into the 
baryon field and the reduced mass by a rescaling, namely 

^b = ,Jl + ^A and rhbia) = mb/{l + ^). Notice 
that rhb is positive definite and only approaches zero as 
a — > oo. 

In the mean-field approximation the Lagrangian gives 



where T3/2 gives the expectation value of the third com- 
ponent of the baryon isospin I^b- The equation of motion 
of the baryons gives the following energy eigenvalue 



Pb = eb = g^b^^ + gpbhbP + \J h^c^kl + 



mld^. (14) 



From the above it is clear that the baryon 6 only exists 
when pb — {gujb^ + gpblzbp) > fhbC^ . The equation of 
motion of the leptons are simply the free Dirac equations 
and are not listed here. Thus, the leptons obey the simple 
relations of a free relativistic Fermi gas: 



ni = ;i7r^, IJ-i=ei = ^h'^c^kf + mfc^. (15) 
For the expectation value of the baryon field we have 

rhbC^k'^dk 



b'i'b 



A 

27r2 



^h^c^k"^ + 
fb rhbC 



27r2 2^3 



hki 



K^kl + TO^c^ 



- 2 2 1 I ^^b 
-rUbC log — 
\ mbC 




and 



{<^b, 



= nb = Jb-^- 



(16) 



(17) 



In the above equations fi = 2Ji -|- 1 is the particle de- 
generacy, while ki and Ci are the Fermi momentum and 
Fermi energy, respectively. 

Since the whole star is considered to be in equilibrium, 
the various processes, such as the direct Urea process 
in Eq. ([2]) and those involving the hyperons in Eq. ([ij, 
impose relations between the chemical potentials of the 
particles. During these processes the baryon number and 
electric charge are conserved, which gives rise to the two 
chemical potentials ps and pQ for baryon number and 
electric charge, respectively. The chemical potential for 
an arbitrary particle with baryon number Bi and charge 
Qi, such as a baryon, lepton or quark, can be written as 
Pi = BiPB+QiPQ- The particle data for all the fermions 
used in this paper are listed in Table. HI 
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TABLE L Particle Data [20]. Here, M is the particle mass, 
B is the baryon number, J is the spin, /a is the 3-component 
ol the isospin, and Q is the charge. 



matter as: 





particle 


M/MeV 


B 


J 


/3 


Q/e 




V 
y 


938.27 


1 


1/2 


1/2 


+1 




n 


939.57 


1 


1/2 


-1/2 







A° 


1115.7 


1 


1/2 










E+ 


1189.4 


1 


1/2 


1 


+1 






1192.6 


1 


1/2 










E~ 


1197.4 


1 


1/2 


— 1 


— 1 


Baryon 


A++ 


1232 


1 


3/2 


+3/2 


+2 




A+ 


1232 


1 


3/2 


+ 1/2 


+ 1 




A° 


1232 


1 


3/2 


-1/2 







A- 


1232 


1 


3/2 


-3/2 


-1 






1315 


1 


1/2 


+ 1/2 









1322 


1 


1/2 


-1/2 


-1 




n- 


1672 


1 


3/2 





-1 


Lepton 


e~ 


0.511 
105.7 






1/2 
1/2 






-1 
-1 




u 


2.4 


1/3 


1/2 


+1/2 


+2/3 


Quark 


d 


4.9 


1/3 


1/2 


-1/2 


-1/3 




s 


105 


1/3 


1/2 





-1/3 



The coupling constants gab, gu>b and gpb for nucleons 
can be fitted with data of the interactions at densities 
near the nuclear saturation point. However, for hypcr- 
ons the values are not known. Here the coupling con- 
stants are simply chosen to be equal among the differ- 
ent baryons [I]. Their values are (5ct/"1ct)^ = 7.487 fm^, 
{gi^lm^f = 2.615 fm^ and {gp/mpf = 4.774 fm^. No- 
tice that the ratio of the coupling constant and the cor- 
responding meson mass is enough to solve this model, 
because we can always rescale the meson fields with their 
masses, so the meson mass will not appear in Eqs. (|10f 

ini). 



Given the chemical potentials and as the in- 
put, together with the neutrality condition ^j^n-iQi = 0, 
Eqs. (|10m3p form a set of self-consistent nonlinear equa- 
tions which can be solved numerically. With the Fermi 
momenta ki of the particles as the output, we can calcu- 
late the energy density and the pressure of the hadronic 




hkb^/h^kf + mlc^{2h^kl + mlc^) 
hkb 



E 



fic 



4 4 1 

m; c log 



hki^h?kf + mfc'^{2h^kf + mfc^) 
hki , 



771/ C 



(18) 



fbC 



hkbJK^kl + mlc:^{2h^kl - Zmlc^) 



\ rribC y mtcr 



1^ 2 2 2 2 2 2 



E 



fic 



hki\lH^kf + mfc^{2h'^kf - 3mfc^) 

(19) 



hk. 



h?kf 

H'- I ^ V ^ ^ 



+3r77fcMog — L + Jl 



The relation between the pressure p and the energy den- 
sity p gives the equation of state, which is shown in Fig.[2l 
and is subsequently used in the TOV equations to get 
the stellar structure. The particle densities as a func- 
tion of baryon density are shown in Fig. |31 Fig. 2] gives 
the particle composition inside the maximum mass star 
with the hadronic equation of state (M = 1.523 Mq, 
R — 9.72 km). Remember that the model is used to cal- 
culate the equation of state for the whole range of the 
density, however, this description may not be appropri- 
ate near the star surface, i.e., the crust. Nevertheless, the 
outer crust only affects low mass stars significantly, while 
we are mainly concerned with maximum mass stars. Fur- 
thermore, in view of the cooling behavior, the outer crust 
has negligible neutrino emissivity compared to the direct 
Urea process in the core. 

Compared to the hadronic phase, the quark phase is 
more easily described. Quark matter is treated as a free 
Fermi gas, whereby we assume that asymptotic freedom 
has taken effect at the very high densities in the center 
of the star. Thus Eq. (|T5|) is also valid for quarks, where 
the lepton chemical potential pi is substituted with the 
quark chemical potential fig and the degeneracy becomes 
fq = 3x {2Jq + 1) = 6 because of the extra color degrees 
of freedom. The fact that the quark phase has a different 
vacuum than the hadronic phase, which has a nonzero 
expectation value of the gluon field, can be taken into 
account by the so-called bag model. This model adds a 
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51E-3 



1200 p 
1000 - 



500 1000 1500 2000 2500 3000^ 



- Hadronic EoS' 
Mixed EoS 

- Nuclear EoS 



p (MeV/fm ) 

FIG. 2: The equation of state of the various matter phases. 
The difference only appears at high density as is shown more 
clearly in the inset, where the dots on the curves indicate the 
central densities and pressures of the corresponding maximum 
mass stars. 




FIG. 3: The particle densities ni/ns versus baryon density 
ub for hadronic matter. 




constant shift called the bag constant Be to the pressure 
and energy density of the quarks, such that 



Pq 



E 



16c7r2/i3 



-m^cMog^ + Jl+ 2 2 



E 



fi 



4 4 1 
-771; C log 



hki^ffkf + mfc^{2h'^kf + mfc^) 
hki 



mic 




Be, (20) 



Pq = 



E 



fq 



' I J 9 



E 



- Be, (21) 




r(km) 



+3771; C log 



where the leptons still contribute because the unequal 
masses of the u, d and s quarks result in unequal densities 
of these flavors, such that even in the pure quark phase 
charge neutrality cannot be satisfied without leptons. 

As discussed in the introduction, there could also be a 
mixed phase between the hadronic phase and the quark 
phase. These two phases are taken to be in equilibrium, 
i.e., the two phases have the same temperature (both set 
to zero here), pressure, and chemical potentials. Charge 
neutrality then determines the volume fraction of these 
two phases. For the three possible phases (hadronic, 
quark, and mixed) the system will be in the one with the 
highest pressure, or equivalently, the lowest grand po- 
tential, which determines the phase transition behavior. 
According to Eq. (PT|) . the bag constant determines the 
scale at which deconfinement sets in. In other words, the 
larger the bag constant, the later the quark phase sets in 
with increasing density. Throughout this paper the bag 
constant is taken to be Be = 230 MeV/fm^ unless indi- 
cated specifically. In Fig. [5] it is shown how the phase 
transition takes place between the two phases. Similar to 
Fig. [3] for the hadronic phase, the various particle densi- 
ties for the mixed phase are shown in Fig. [S]as a function 
of baryon density. Fig. [7] shows the particle composition 
inside the maximum mass star with the mixed equation 
of state (M = 1.479 Mg, R = 9.81 km). 



C. Thermal evolution equations 



FIG. 4: The particle composition of the maximum mass neu- 
tron star with the hadronic equation of state. 



The equations governing the thermal evolution of a 
spherically symmetric star, given by the metric in Eq. 
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FIG. 5: The pressure surface of the hadronic phase (dark 
gray) and the quark phase (hght gray). The thick sohd hues 
are the neutrahty curves in each phase. The thick dashed 
curve is the intersection of the two pressure surfaces and shows 
when the two phases are in equilibrium. The arrows indicate 
how p and hq change with increasing /ib ■ 




0.01 



FIG. 6: The particle densities Ui/nB versus baryon density 
ns for the mixed phase case. The densities of the quarks are 
expressed in terms of their baryon densities instead of their 
number densities. 



are ml: 



cv 



dt 
dr 



-X 



47rr^ 



9?" 



2* 



-,2<S> 



(22) 
(23) 



where cy is the specific heat capacity at constant volume, 
K is the thermal conductivity, and q-y are the neutrino 
and photon emissivity, respectively. The photon emissiv- 
ity q-y is only nonzero on the surface r = R. It is con- 
venient to define the redshifted temperature T — Te* 
inside the star. Similarly, Le^* represents the redshifted 
luminosity corresponding to the heat current. In order to 



FIG. 7: The particle baryon densities inside the maximum 
mass neutron star with the mixed equation of state. Note 
that the quarks have baryon number 1/3. 



compare with observations, the photon luminosity is of 
great importance. Here it is assumed to obey the black 
body radiation law Lj = AnR^q^ = AiraR^Tf, where 
(J = 5.67 X 10^^ erg/cm^K^s is the Stefan-Boltzmann 
constant and Tg is the surface temperature. According 
to general relativity, an observer at infinity will measure 
the gravitationally red-shifted temperature and luminos- 
ity. The observational quantities are thus Too = TsC**-^^ 
and Loo = L-ye^**-^^. As is generally accepted, a very 
thin layer of the outer crust of the neutron star will act 
as a thermal insulator causing the temperature at the sur- 
face to be much lower than inside the star. The relation 
between the surface temperature and the inner tempera- 
ture depends on the chemical composition of this envelop. 
We will simply locate such a layer at the surface, neglect- 
ing its thickness, and adopt the T^o-T relation given by 
Potekhin et al. [24I, as shown in Fig. |51 and set T to be 
the inside temperature at r = i?. Note that the surface 
temperature Ts introduced above is the outside temper- 
ature bX r = R. 

The parameters cy, k, q^ and q^ need to be determined 
to solve the thermal evolution equations, Eqs. ([22]) and 
(|23l) . The specific heat cy can simply be obtained by 
summing up all the contributions of the different fermions 
in the model. In the low-temperature approximation, 
Fermi-liquid theory gives 



CVi 



(24) 



IS 



where fcs is Boltzmann's constant, m* = hki/v^ 
the effective mass on the Fermi surface, and Vi = 
dei/ (hdk)\k=ki is the Fermi velocity of quasi-particle i. 
For leptons the effective mass is easily obtained since 
they are described as a free Fermi gas, i.e., m^* ~ 
y^h?kf /c^ + mf = fii/c^. Similarly for quarks, m* = 
Hq/c^- For baryons the effective mass follows from the 
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FIG. 8: The T^o-T relation for different numbers of ligiit el- 
ements parameterized hy rj = gi4^AMi/M, where AMi is the 
mass of the light elements in the envelop and (714 is the sur- 
face gravity in unit of 10^'* cm/s^ [S^l, the different curves 
correspond to log?; = —00, —16, —14, —12, —10, —8, —6 and 
from bottom to top, respectively. These results are for the 
maximum mass neutron star with the hadronic equation of 
state, where gi4 = 2.918. 



dispersion relation in Eq. ([141 



Vb 



such that mt 



deb{k) 



hdk 



hkh 



k—ku 



ml' 



(25) 



/c^ + ml. 



The effect of interac- 



tions is, in the mean-field approximation, only present in 
the reduced mass rhb of the baryons. 

The general expression for the neutrino emissivity has 
the form q^i = CiT", where Ci and s are different con- 
stants for each kind of process. For the direct Urea pro- 
cess we have s = 6 and for the modified Urea process 
s ~ 8. The difference in the exponents shows the ineffi- 
ciency of the modified Urea process since the temperature 
T of a neutron star is much smaller than the typical Fermi 
temperature Tp. The neutrino emissivities for the vari- 
ous processes are summarized in Ref. Q . For the general 
baryon direct Urea process, cf. Eq. the emissivity is 



q^i2i « 1.207 X 10^ 



(1 MeV)TO2 



R12TIQ121 erg/cm^s, 



(26) 

where Tg is the temperature in units of 10^ K, Q121 is 1 
only if the Fermi momenta of the two baryons hi , &2 and 
the lepton I can form a triangle, otherwise 612; — 0. The 
coefficients R12 vary depending on the baryons involved 
in the process. In our mean-field model the maximum 
mass star contains massive hyperons up to S, such that 
all possible direct Urea processes given in Ref. [2^ need 
to be included, as listed in Table HIl For the quark direct 



TABLE II: The coefficients R12 for different direct Urea pro- 
cesses [2^, denoted by the baryons involved in each process. 



np 


A«p 


E"n 


E-A° 




H-A" 


H-E" 


H°E+ 




1 


0.0394 


0.0125 


0.2055 


0.6052 


0.0175 


0.0282 


0.0564 


0.2218 



Urea emissivities, we refer to Ref. |24| : 



qvudi ~ 4.773 X 10 
q^usi ~ 2.552 X 10 



1 McV^ 



Tl erg/cm^s. 



l^h^c^[{ku + klf-k1]^l, . 



1 MeV^ 



erg/cm'^s, 
(27) 



where the efficiency of the strange quark channel is re- 
duced because it changes strangeness. Since the direct 
Urea channels are open, it is not necessary to consider 
other neutrino-emission processes. 

In the present model, the thermal conductivity k can- 
not be easily obtained since it strongly depends on the 
interactions between the particles in the complicated 
hadronic phase. For the mixed case, the interface be- 
tween the phases may introduce even more uncertainty. 
Therefore, the thermal relaxation period inside the star 
is not considered here and the calculation simply starts 
when the star has become isothermal, as indicated in the 
introduction. Because of this simplification the inter- 
esting early stage behavior of the cooling cannot be dis- 
cussed. This disadvantage may also be one of the reasons 
why such hadronic models are not extensively studied. 
However, the typical thermal relaxation time is usually 
less than a century and most of the available observa- 
tional data is of the later stages, i.e., after the star has 
become isothermal, such that the results obtained after 
the relaxation period has ended are still useful. Conse- 
quently, the star will be treated as an isothermal object 
with a constant temperature T and the equations are 
rather simplified. Instead of two coupled partial differ- 
ential equations in Eqs. and we have now only 
one ordinary differential equation: 



at 



(28) 



where the capital letters represent the integrated param- 
eters over the volume of the star, namely 



,(s-2)*(r) 



(29) 
dr, (30) 



where the T-dependent parts are taken outside of the 
radial integral and the temperature-independent prefac- 
tors of cvi and qi,i are denoted with a bar. The factors 
of e*^'"-' in the denominators are a consequence of the 
construction of the isothermal temperature T. Note that 
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the isothermal temperature does not take into account 
the temperature decrease near the surface of the star, 
which occurs in a very thin layer close to the surface. 

The total heat capacity and the neutrino emissivities 
for both the direct and modified Urea processes as a func- 
tion of the radial coordinate are shown in Figs. |9] and 
1101 for the case of a maximum mass neutron star with 
the hadronic and mixed equations of state, respectively. 
Here "total" means summed over all particles which con- 
tribute, but not integrated over the star volume. The 
quantities were obtained at T = 10^ K, but the values 
at different temperatures can be obtained directly from 
the temperature dependence of cv and q^. The discon- 
tinuities in the emissivity are due to the step functions 
in the expressions for the emissivities. These step func- 
tions are a consequence of the low-temperature approx- 
imation, which docs not take into account any further 
momentum fluctuations around the Fermi surfaces. Also 
^MUrca IS Calculated only for the nucleon processes shown 
in Eq. ^ with the following emissivity [3|: 



1.882 X 10 



(luMn 



quMp ~ — - 
V ml 



19 



(1 MeV) fiimf^mp 



erg/cm^s. 



-\- hi kj^ 
8k„ki 



-QupqvMm (31) 



where the subscript Mn or Mp means the process is mod- 
ified by a bystander neutron n or proton p, and is 
1 only if Zkp + ki > kn as required by momentum conser- 
vation in the Mp process. The modified Urea processes 
are only calculated when the nucleon direct Urea pro- 
cess is forbidden, since the formulae are only valid in this 
case. The small region of overlap where both gourca and 
9MUrca are nonzero is due to the different thresholds for 
the direct Urea processes with electrons e~ and muons 



"^DUrca 




r(km) 



FIG. 9: The total heat capacity (cv) and the neutrino emis- 
sivities of the direct Urea (gourca) and modified Urea (gMUrca) 
processes inside the neutron star with the hadronic equation 
of state at T = 10® K. Here, cy is in units of 10^^ erg/cm'^K, 
QDUrca IS in units of 10^^ erg/cm^K, and giviurca is in units of 
10^® erg/cm^K. 




III. NUMERICAL RESULT 

In this section the three different equations of state for 
the hadronic matter, the mixed phase of hadronic and 
quark matter, and nuclear matter are compared. Since 
the maximum masses are not very large for the last two 
cases, we concentrate on the star with the maximum al- 
lowable mass for each equation of state. The initial tem- 
perature is always set as T = 10^ K at t = 0, which 
corresponds to the time at which the isothermal condi- 
tion is reached. As mentioned previously, the thermal re- 
laxation period is not included in our calculation, which 
is about several decades according to calculations with 
nuclear matter 0]. The results are plotted in a loga- 
rithmic time scale, such that our results can be shown 
directly with most of the observational data even though 
the data include the thermal relaxation era of the neutron 
stars, that is negligible compared to the age of the stars. 
However, this makes it hard to draw a comparison with 
the recent observation of the cooling of the Cassiopeia A 



FIG. 10: The same as in Fig.[9]but now for the mixed equation 
of state. 



neutron star [25|, which provides direct evidence for the 
cooling process. This neutron star is too young and is 
believed to have become isothermal quite recently, such 
that the thermal relaxation period is not negligible. 

In order to avoid confusion, our results always show 
cooling curves starting at t = 1 yr, when the result has 
become not very sensitive to the initial conditions. In 
fact, the effect of the initial temperature is limited to the 
very beginning. The higher the initial temperature is, 
the faster the dependence dies away, which can be seen 
clearly from the asymptotic temperature dependence of 
the neutrino-emission dominated era 



is-2)Qc 
Co 



(32) 



where Co and Qq are the constant coefficients of the T- 
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dependent factors in Eqs. ([29]) and (pO)) . This asymptotic 
solution can be easily obtained by dropping the photon 
radiation term in Eq. ([28|) and keeping only the leading 
neutrino-emission process with the smallest exponent s. 
Note that the exponent s is never smaller than 6 for any of 
the neutrino-emission processes. It is quite clear that the 
late-time behavior of the neutrino-emission dominated 
era is thus completely determined by the coefficients of 
the heat capacity and the neutrino emissivity and not by 
the initial condition T{0). 

For the hadronic equation of state, Fig. [TT] shows the 
evolution of the temperatures T and Toe, where it can 
be seen that the cooling process can be divided into 
two stages, namely the neutrino-emission dominated and 
photo-radiation dominated era. The two eras can be 
clearly seen from the energy loss due to the different pro- 
cesses as shown in Fig. ll2[ where the switch from neutrino 
emission to photon radiation occurs at t « 3 x 10^ yr. In 
Fig.[T2]the energy emission due to the modified Urea pro- 
cesses is also shown demonstratively by considering only 
the processes involving the nucleons in Eq. ([3]). This, of 
course, underestimates the energy loss due to the modi- 
fied Urea processes, however, the magnitude is expected 
to be of the same order. The energy loss by these pro- 
cesses is seen to be less than 10~^ of the total energy 
loss, such that the modified Urea processes are always 
negligible in the present calculation. In Figs. [TT] and [T^] 
the Too-T relation is used with 7] = 10^^". It is found 
that different values of rj only slightly change the cooling 
process, as shown in Fig. 1131 where two groups of cool- 
ing curves with extremely large and almost vanishing ry 
are compared. The difference is that, for the neutrino- 
emission dominated era, Too can differ by about a factor 
of 2 but the inner temperature is almost not infiuenced. 
With less massive elements in the envelop, i.e., smaller 77, 
the turning point into a photon-radiation dominated era 
is a little sharper and the final temperatures are a little 
lower. In any case, these effects are not very significant 
and can hardly be distinguished with the present accu- 
racy of observational data. Without loss of generality, 
the value 77 = 10"-'^° will be used in this section. Here 
we should also point out that the calculation cannot be 
carried out after T is smaller than 10"* K. because then 
the Too-T relation given in Ref. |22] is no longer valid. 

The thermal evolution of the neutron star with the 
mixed equation of state is shown in Fig. 1141 which is 
similar to the purely hadronic case. The straight lines 
represent the asymptotic behavior, which are solved us- 
ing Eq. by neglecting or for the neutrino- 
emission dominated or photon-radiation dominated era, 
respectively. For example, in the neutrino-emission dom- 
inated era, for the hadronic equation of state the tem- 
peratures scale with time as T ^ 1.6 x 10^ (t/yr)~^/^K, 
Too ^ 1.1 X 10'^ (t/yr)^0-602/4K, where the power 0.602 
comes from the asymptotic approximation for the Too- 
T relation at low temperature, Too oc For the 

mixed phase, the temperature scalings are just slightly 
higher: f - 1.77 x lO^ {t/yry^/'^K, Too - 1-13 x 




Log t (yr) 

FIG. 11: The temperature evolution of the neutron star with 
the hadronic equation of state. The solid curve represents 
T while the dashed curve represents Too, with the dotted 
straight lines indicating the asymptotic behavior. 
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FIG. 12: The energy lost by different processes in the neutron 
star with the hadronic equation of state, where the solid curve 
is the total energy loss rate, the dashed curve is the energy 
loss due to direct Urea neutrino emission, the dotted curve is 
due to modified Urea neutrino emission, and the dot-dashed 
curve is due to photon radiation. 



10^ (t/yr)-°-6°2/4j^^ rpj^jg .g ^^^^ ^j^^^ ^j^^ 

mal parameters, after being integrated over the whole 
star, are not quite different for the two cases even though 
the stars have different structure and mass. For ex- 
ample, for the hadronic case we have Cy = 3.331 x 
W^Tq erg/K and = 1.227 x 10'^'' (fg)^ erg/s, while 
for the mixed phase Cy = 2.839 x lO^^Tg erg/K and 
= 7.280 X 10'*'^(T'9)^ erg/s. To demonstrate the cflFect 
of the mixed phase, the same curves are shown in Fig. 1151 
with a smaller bag constant of Be = 170 MeV/fm'^, where 
we see that it causes a slightly higher temperature. As 
pointed out earlier, a smaller bag constant causes the 
deconfined phase of quarks to appear at lower densities, 
such that the volume of the mixed phase is increased. 
In this case we find Cy = 2.156 x lO^^Tg erg/K and 
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Log t (yr) 

FIG. 13: The temperature evolution of the neutron star with 
the hadronic equation of state with different Too-T relation. 
The solid curves correspond to r; = 1 while the dashed curves 
correspond to r; = 10"^". The upper group is T and the lower 
group is Too. 




Log t (yr) 

FIG. 14: The same as in Fig. Illl but with the mixed equation 
of state. 



= 2.590 X 10'*^(f9)6 erg/s. Nevertheless, the total 
effect on the cooling process is not drastically changed. 
According to these results, we expect that even if the 
mixed phase is reduced because of screening and surface 
effects, the cooling behavior is not changed considerably. 

For the nuclear equation of state, the maximum mass is 
1.719 Mq with R = 10.04 km and its thermal evolution 
is shown in Fig. 1161 To summarize, the luminosity for the 
three types of equations of state arc plotted in Fig. [T71 
where a comparison with observational data is also made. 
In all three cases the cooling is too fast compared with 
the data. The star with the nuclear equation of state 
cools a little faster, since the direct Urea processes with 
hyperons in the hadronic phase and the mixed phase are 
not as efficient as the nucleon direct Urea process [23| . 
However, the difference is quite small such that it might 
be hard to distinguish these different equations of state 
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FIG. 15: The same as in Fig.[Tlbut with = 170MeV/fm^ 

10 % 
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FIG. 16: The same as in Fig.[TT]but with the nuclear equation 
of state. 



simply by their thermal evolution after the isothermal 
condition is reached. 



IV. EFFECT OF SUPERFLUIDITY 

From the above results, we can see that our neutron 
stars cool too fast because the direct Urea process is open 
for all the three types of equations of state. An expla- 
nation could be that the neutrino emissivity is overesti- 
mated by neglecting the possibility of superfluidity. As 
pointed out in the introduction, it is generally believed 
that superfluidity appears inside neutron stars due to the 
attractive part of the nuclear force, for both the nucle- 
ons and the hyperons. Due to the presence of an energy 
gap A, associated with the binding energy of the Cooper 
pairs, the number of excitations near the Fermi surface 
are suppressed by a factor of exp(— A/fc^T) when the 
temperature is smaller than a certain critical tempera- 
ture Tc- Therefore, pairing reduces the heat capacity 
and the neutrino emissivity signiflcantly when the star 
cools down to temperatures below thus the cooling 
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FIG. 17: The luminosity of neutron stars with different equa- 
tions of state, with solid, dashed, and dotted curves repre- 
senting the hadronic, the mixed, and the nuclear equation of 
state, respectively. The observational data are from Ref. [2^ . 
with the numbers indicating the corresponding stars as: 1 - 
RX J0822-4247, 2 - IE 1207.4-5209, 3 - RX J0002-)-6246, 4 - 
PSR 0833-45 (Vela), 5 - PSR 1706-44, 6 - PSR 0538-^2817, 
7 - PSR 0656-hl4, 8 - PSR 0633-H748 (Geminga), 9 - PSR 
1055-52, 10 - RX J1856. 5-3754, and 11 - RX J0720.4-3125. 

behavior is drastically changed. In this section we will 
explore the effect of superfluidity to the cooling process 
of neutron stars. 



A. Superfluidity inside neutron stars 

In our present model for neutron stars, there are dif- 
ferent fermions, namely the baryons, the leptons, and the 
quarks. The typical Tc of baryons is about 0.1 — 1 MeV, 
which is close to the initial temperature of neutron stars, 
such that the superfluidity of baryons is very impor- 
tant to the cooling process. As for the leptons, such 
as electrons, they could become superfluid due to in- 
teractions via phonons, however, the typical Tc is very 
small, namely in the order of several kclvin. Hence lep- 
ton superfluidity is expected to be unimportant to the 
neutron star cooling process. For quarks in the mixed 
phase, the gaps can be about 50 — 100 MeV [13], which 
is much higher than the typical temperatures of the stars 
we are concerned with. The quark contribution will thus 
be suppressed more strongly than the baryonic contri- 
bution. Without superfluidity the quark contributions 
to the heat capacity and the neutrino emissivity were 
shown to be of the same order as the baryonic ones, such 
that with superfluidity they are completely negligible at 
the same stellar temperature. Throughout this section 
we thus only consider the contributions to cy and of 
non-superfluid leptons and superfluid baryons but neglect 
those of the superfluid quarks. 

Superfluidity not only reduces the neutrino emissivity 
but also opens a different channel of neutrino emission 
based on the breaking and formation of Cooper pairs. At 



temperatures not far below Tc thermal fluctuations can 
cause pair break-up into single-particle excitations, which 
subsequently reform into pairs. Neutrino emission due to 
such processes is quite efficient at temperatures slightly 
below Tc and can even surpass the direct Urea process in 
some cases [23|. However, its contribution also decreases 
dramatically when T becomes small due to the exponen- 
tial reduction factor mentioned previously. According to 
a similar argument as above, we only need to consider 
such neutrino-emission processes involving baryon pairs. 

By taking into account superfluidity, the dominant 
neutrino-emission process is no longer simply determined 
by the density, which is crucial for the threshold of the 
direct Urea process, but also depends on the reduced tem- 
perature T = T/Tc- However, it is still true that the mod- 
ified Urea process can be neglected as long as the direct 
Urea process is open, since the modified Urea process in- 
volves more particles and is more strongly suppressed due 
to the energy gaps of all superfluid participants. Thus we 
again do not consider the modified Urea process since the 
direct Urea process is always open in our study. Accord- 
ing to Ref. [27[ electron-electron bremsstrahlung becomes 
a dominant process when all baryons are strongly super- 
fluid at low temperature. So for a correct estimate of the 
neutrino emission at low temperatures we also include 
lepton-bremsstrahlung processes. Based on the discus- 
sion above, in order to include the effect of superfluidity 
we should recalculate cy and for all baryons, neglect 
those of the quarks in the mixed phase, keep the lepton 
contributions unchanged, and also include the newly in- 
troduced neutrino emission due to baryon pair break-up 
and lepton bremsstrahlung. We next discuss these effects 
separately. 

To describe the superfluidity of the particles, it is nec- 
essary to specify the type of pairing. At low density, 
the singlet-state nuclear interaction is attractive, while 
at higher densities it becomes repulsive. It is believed, 
however, that the triplet-state interaction still provides 
an attractive channel at high densities, such that triplet- 
state pairing is expected in the core of the neutron star. 
This transition happens around the nuclear saturation 
density no = 0.16 fm"'^. The proton pairing is usually 
taken to be in the singlet channel, even in the stellar 
core, due to their low concentration. However, in our 
model, as seen from Figs. H] and [3 the densities of neu- 
trons and protons in the core are quite close to each other, 
namely rip ^ 0.2 — 0.3 fm~'^ and n„ ^ 0.4 fm~'^. It thus 
seems that both the protons and the neutrons should 
form triplet-state pairs. The hyperons, which can also 
become superfluid, are usually taken to pair in the sin- 
glet channel. But in our model there are some hyperons, 
e.g. A", which can have densities comparable to the neu- 
trons and protons. As was discussed earlier, the details 
of the interactions between the hyperons are not well es- 
tablished and thus there are some ambiguities in dealing 
with hyperons with largely varying densities. For the 
sing let-state pairing of neutrons, protons [2^, and A°'s 
[28j the results of many models largely agree on the order 
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of magnitude of the pairing gaps and their density depen- 
dence. As discussed, triplet-state pairing is of great in- 
terest to our model, however, there is little known about 
this kind of pairing except for the case of neutrons for 
which it is still highly model dependent. In fact, the ex- 
istence of triplet-state pairing inside neutron stars is still 
uncertain. According to observations of the cooling neu- 
tron stars, it seems necessary to have all baryons in the 
superfluid state. For example, the effect of superfluidity 
on the thermal evolution of neutron stars with hyperons 
in the core has been studied by Schaab et al. [2^ , which 
only included the singlet-state pairing of A" at low densi- 
ties but did not include the triplet-state pairing of A" at 
high densities and ignored the pairing of other hyperons. 
The cooling was found to be too fast for heavy stars since 
not all the direct Urea processes were suppressed. 

To estimate the effect of superfluidity, we make the 
simplifying assumption that, regardless of the density, 
the neutrons, protons, and A'^'s pair in the triplet chan- 
nel, the remaining hyperons pair in the singlet channel, 
and the critical temperatures of all the baryons are taken 
to be equal. Equating all critical temperatures of the 
various baryons is consistent with our model, since all 
baryons couple equally to the meson fields. Although our 
pairing mechanism seems quite crude, we expect that, at 
least qualitatively, the effect of pairing on the thermal 
evolution of the star will be taken into account. To dis- 
cuss our simplifications, we first need to go into the de- 
tails of how the specific heat and neutrino emissivity are 
reduced by superfluidity. 



1. Reduction factor 

The new expressions for cy or are easily obtained by 
multiplying the original expressions by a reduction factor 
Rc or Rq. These reduction factors are functions of the 
reduced temperature r and also depend on the type of 
superfluidity considered. A rather systematic calculation 
of these reduction factors has already been carried out 
(see review (2^ and references therein). In general, the 
different types of triplet-state pairing can be represented 
by the projections mj = 0, ±1,±2 of the total angular 
momentum. We only present the most studied mj = 
case, which is denoted as type-B pairing. Singlet-state 
pairing is referred to as type A. By introducing the di- 
mensionless variables r = T/Tc and v{t) = A{T)/kBT, 
the properties of superfluidity can be described indepen- 
dent of Tc- Numerical fits for the energy gap are given 

by 



.. = 71^(1.456-^ + 1:1^), (33) 
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FIG. 18: The reduction factors Rc (solid curves) and Rq 
(dashed curves) for different types of superfluidity, where 
type-A (thick curves) is always lower than type-B (thin 
curves) pairing. Notice the jump of Rc a.t t = 1. The dotted 
curve represents R^g used for the direct Urea process where 
both baryons are superfluid and have the same Tc- Also RqAA 
[30t (dot-dashed curve) and i?^^ (thick dot-dashed curve) are 
shown for comparison. 



and u = at r > 1. The reduction factors for cy are 

, 2.5 

RcA = ( 0.4186 + Jl.014 + 0.2510^2^ ' 



RrR — 



^ gl. 456-^2.120-1-1)^ 

0.6893+ a/0.6241 
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X e 



1.934-y^3?740T^ 



(35) 



(36) 



Notice that Rca = 2.426 and Rcb = 2.188 are greater 
than 1 at T = 1, such that the specific heat is discon- 
tinuous as the temperature falls below Tc- The neutrino 
emissivity for the direct Urea process is different if both 
or only one of the two involved baryons is superfluid. For 
the case with only one superfluid baryon 
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Generally, type-A pairing always has a stronger reduction 
effect than type B, as shown in Fig. [THl 

The reduction factor for the case where both particles 
are superfluid is not simply the product of the two reduc- 
tion factors of each particle, due to different phase-space 
restrictions. This was pointed out in Ref. (soj . where 
the authors presented the numerical fits for the reduc- 
tion factors for the AA and BA cases. Here the cases are 
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distinguished by the type of pairing (A or B) of the neu- 
trons and protons, labeled as AA, BA or BB, where the 
first letter signifies the type of pairing for the neutrons 
and the second the pairing type for the protons. We do 
not present their complicated expressions here. The most 
important factor we need is for the BB case, which as far 
as we know has not been explicitly calculated yet. There- 
fore, we use RqB{Ti)RqB{T2) instead of RqBB{Ti,T2) for 
those direct Urea processes with both baryons paired in 
the triplet-state, which we refer to as the standard set- 
ting. Furthermore, we even use this expression for all 
the direct Urea processes with two superfiuid baryons, 
independent of the pairing type. Therefore, the effect of 
superfiuidity on the emissivity for the BB-type process is 
expected to be overestimated [l^]. On the other hand, 
since the above substitution is also applied to the n, p 
and A*^ at lower densities and other hyperons, which are 
supposed to pair in the singlet channel, the suppression 
will be underestimated in these cases. The errors in these 
two approximations thus work in opposite direction and 
roughly reduce the total error of our calculation. The re- 
duction factors of the direct Urea process with two super- 
fiuid baryons are compared in Fig. 1181 It is seen that their 
differences can become very large for T <^ Tc- However, 
in this case the reduction factors are extremely small and 
the lepton-bremsstrahlung processes are expected to be 
dominant, such that in this regime any ambiguity due to 
the type of pairing is expected to be unimportant. 



TABLE III: The factor aj, of the neutrino emissivity due to 
the Cooper pair break-up of the various particles [Sll]. Notice 
that we take the coefficient of A'' to be the same as n, which 
differs from Refs. 
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Of, for the various baryons used in the calculation is shown 
in Table Uni 

The function F{v) plays the same role as the reduction 
factor and is fitted in Ref. [311 to: 



Fa =(0.602wi + 0.5942< + 0.288O ( 0.5547-h 

1/2 



0.1983 + 0.0113u^ 



„2. 245-^5.04-1-41)2 



(40) 



Fp, 



1.204W?, + 3.733i;i 



0.3191w| 



0.3511w| 



0.7591- 



0.05803 0.3145w| 



^0. 4616-^0. 2131+4i;| 



(41) 



Furthermore, the neutrino emissivity due to the lepton- 
bremsstrahlung process ee is given by [ssj 



2. New neutrino emission processes due to superfluidity 

The neutrino emissivity associated with Cooper pair 
break-up is given by [IH 



1.17 X lO^^iV^ 



mt,c 



abF{v)Tl erg/cm^s, (39) 



where Ni, = 3 is the number of neutrino fiavors, ah is 
a numerical factor from the electroweak neutral currents 
and depends on the quark composition of the baryon and 
the pairing type. As far as we know, no calculation has 
been carried out which included the triplet-state pairing 
of hP . We simply use oa = a„ , considering that the quark 
composition of A" is similar to the neutron with one d 
quark replaced by one s quark, while their contributions 
to the neutral current are the same. This treatment is 
different from the approximation in Ref. [11], where the 
contribution from quarks other than u and d is neglected. 
This may cause some uncertainties, but what matters in 
the cooling process is the order of magnitude of each 
process. Since af, appears as a multiplier rather than 
an exponent in Eq. p9|) . such an inaccuracy will not be 
magnified during the calculation. In fact, we will see 
that due to the assumption of a uniform Tc, the neutrino 
emission associated with Cooper pair break-up is always 
negligible compared to the direct Urea process even if the 
suppression due to superfluidity is included. The factor 



C = 2.089 X 10" 



hck„ 



IMeVy, 



■Tg erg/cm^s. 



(42) 



where ys is a dimensionless parameter representing the 
effect of screening in the plasma, ys = kgc/'^.ke, with ksc 
the screening wavenumber. This wavenumber is obtained 
for the case of static screening in the limit of zero tem- 
perature (3^ by the Thomas-Fermi expression: 



=2 



4e 

I b 



(43) 



where the summation over b is over all charged baryons, 
and Zb represents the effect of superfluidity on the 
baryons. Since lepton bremsstrahlung only becomes 
important if all the baryons are highly superfluid, for 
which Zb becomes negligibly small compared to the lep- 
ton terms, we can simply omit the terms related to the 
baryons. Subsequently, ys is simplified to 



irhckl 



(44) 



where a is the fine-structure constant. Other pro- 
cesses involving non-relativistic muons are calculated in 
Ref. (ssj . In our model the muons are relativistic, so we 
adopt the similarity criterion Q and get 



'iV ^2 



(45) 
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FIG. 19: The neutrino emissivities of the various processes at 
different temperatures. The solid curve is for tire direct Urea 
process, the dashed curve for the triplet-pair break-up pro- 
cess, the dotted curve is for the singlet-pair break-up process, 
while the dot-dashed curve for the lepton-bremsstrahlung pro- 
cess. The vertical dashed line at t = 1 indicates the onset of 
superfluidity of all baryons. The emissivities are calculated at 
the central density of the maximum mass neutron star with 
the hadronic equation of state and Tc = 10^ K. 

To summarize, in Fig. [19] we present the neutrino emis- 
sivities of the discussed processes as a function of reduced 
temperature r after the inclusion of superfluidity. Due 
to our assumption of equal Tc, the emissivity associated 
with Cooper pair break-up is always about 5 orders of 
magnitude smaller than for the direct Urea process. This 
agrees qualitatively with the result in Ref. [l^] for the 
equal Tc case of nuclear matter. For the neutron star with 
the mixed equation of state, the situation is similar since 
the leading contribution of triplet-state paired baryons 
changes very little. In fact, we can always neglect the 
contribution to the neutrino emissivity associated with 
Cooper pair break-up. 

Finally, it should be pointed out that superfluidity can 
also affect the equation of state. However, the equation 
of state used in the TOV equations is an integral over 
all states in momentum space and it is expected that the 
effect of the energy gap at the Fermi surface is very small 
for very high densities. Therefore, we can still use the 
stellar structure obtained from the previous unaffected 
equation of state. Again, we only present the results 
from the maximum mass star for each type of equation 
of state. 



B. Numerical results 

The calculation of the thermal evolution of the star 
is similar to the case without superfluidity, but now the 
coefficients Cy and Qi, in Eqs. (P^)) and include the 
temperature-dependent reduction factors. Note that the 
relevant parameter for the thermal evolution is T rather 



than T, such that the reduction factors are functions of 
r(r) = T'e~*(''V7c, which is not constant throughout 
the star if we use a constant Tc. This radial dependence 
should be included in the integration of Eqs. and 
(pO| . Since the r dependence of the reduction factor is 
nonlinear, the variables T and r cannot be separated, 
which means that the numerical integration must be car- 
ried out at each instance to solve the differential equa- 
tion, cf. Eq. (j28p . To circumvent this we introduce an 
additional approximation by taking Tee*'-'"-* a constant, 
such that the reduced temperature r is constant inside 
the star. This decouples r and r in Eqs. (PH]) and (|30p . 
which means we do the radial integration only once. 
Then Eq. ([25)1 is independent of the radial coordinate 
and the reduction effect is represented by several extra 
T-dependent terms in the coefficients. As is well known, 
the critical temperature Tc is a function of the density, 
however, the dependence is often uncertain. Because of 
the discrepancies between the pairing models, it is even 
unclear whether Tc for each type of baryon increases or 
decreases with baryon density in the range of interest. 
In fact, setting Tc = const throughout the star also im- 
plies some kind of density dependence because the baryon 
density changes with the radius. In other words, as an 
ansatz. Tee*'*"' — const is as reasonable as Tc = const, 
even though the former seems unnatural because of its 
dependence on the macroscopic stellar structure. We use 
the former ansatz first because of its simplicity. Later 
we make a comparison between these two ansatze. In 
any case, e"*^''-' only varies smoothly and monotonically 
around 2, for example, from about 2.2 at r = to about 
1.4 at r = i? for the maximum mass neutron star with 
the hadronic equation of state, so we expect the difference 
between these two ansatze will not be of several orders 
in magnitude. 

In the numerical calculation with the first ansatz, we 
set Tce*('') = 'f(O) = 10^ K such that the baryons are 
supcrfluid from the start. This implies we have Tc ~ 
2.2 X 10^ K at the central density while Tc ~ 1.4 x 10^ K 
near the surface. If, however, r(0) were to be set higher, 
the reduction factors would at first not take effect and the 
direct Urea process would thus quickly reduce T down to 
Tc after which it would be suppressed due to the onset of 
superfluidity and the cooling slowed down. This process 
takes only several seconds. On the other hand, if T(0) 
would have been slightly lower, for example, above O.lTc, 
then in the beginning the superfluid suppression would 
not be very strong and the direct Urea process could 
still cool the star to low temperature, i.e., the highly su- 
perfluid regime, within one year. Since the isothermal 
assumption we used here is virtually a later stage behav- 
ior, given a reasonable Tc and the well-accepted range 
of T{0) of neutron stars, the result we presented in this 
section depends little on the accurate values of the tem- 
perature. As before, we take the neutron star with the 
hadronic equation of state as the basic example, then we 
study the stars with the mixed and nuclear equation of 
state. 
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FIG. 20: The luminosity of a neutron star with the hadronic 
equation of state. The solid, dot-dashed, dashed, and dotted 
curves are for 77 = 10~®, 10~^^, 10~^^, and 10~^^, respectively. 
The observational data are again shown for comparison. 



We find several general properties of the cooling be- 
havior of stars containing superfluid matter from the nu- 
merical results. One obvious difference with the non- 
superfluid case is the shape of the cooling curve. Since 
the neutrino-emission processes, which dominate the 
early stage cooling, are strongly suppressed with decreas- 
ing temperature, the temperature stays high for a longer 
time and the shift to the photon radiation era takes place 
at higher temperatures. In contrast to the large slope 
in the non-superfluid case, we find a rather slowly de- 
creasing plateau around T ^ 10^ K and Too ^ 10^ K, 
followed by a much sharper turn to the photon radia- 
tion stage as can be seen in the following figures. An- 
other difference can be seen in that the parameter rj in 
the Too-T relation plays a more significant role. Because 
the neutrino-emission processes are suppressed due to su- 
perfluidity, photon radiation becomes more important at 
higher temperatures and the dependence on how the in- 
ner temperature is screened by the surface layer becomes 
more visible. We present the luminosity curves with dif- 
ferent 77 in Fig. [201 where we see that the observational 
data favor moderate values of 77. In the following, we 
usually use 77 = 10~i^"^. 

Here some other general properties are summarized. 
First, we find that the heat capacity cy of the baryons is 
quickly reduced as a consequence of superfluidity. Start- 
ing from t w 1 yr, the heat capacity is completely dom- 
inated by the Icpton contribution, such that the pair- 
ing type of the baryons is seen to be of no importance 
to the heat capacity. Second, as was already shown in 
Fig. 1191 neutrino emission via Cooper pair break-up is 
always negligible compared to the direct Urea process. 
We also find that the role of lepton bremsstrahlung in 
the cooling process is negligible since in the later stages 
when it dominates the neutrino emission, photon radi- 
ation has already become the dominant cooling effect. 
Therefore, the cooling process is determined only by the 



FIG. 21: The energy lost by the various processes in the neu- 
tron star with the hadronic equation of state, where the thick 
solid curve is the total energy loss rate, the thin solid curve 
is for the direct Urea process, the two dot-dashed curves are 
for the Cooper pair break-up processes with the upper one 
for the triplet-state pairing and the lower one for the singlet- 
state pairing, the dotted curve is for lepton bremsstrahlung, 
and the dashed curve is for photon radiation. 




Log t (yr) 

FIG. 22: The temperature evolution of the neutron star with 
the hadronic equation of state. The various curves represent 
different pairing types. The solid curves are for our standard 
setting, the dotted curves are when all baryons pair in the 
singlet-state, and the dashed curves are also for singlet-state 
pairing but with the reduction factor Raa replaced by i?^. 
As before, the upper group is T and the lower group is Too, 
respectively. 



direct Urea process and photon radiation. The various 
energy loss rates arc shown in Fig. 1211 Third, the direct 
Urea process depends on the type of pairing. However, 
this dependence is quite moderate. We show the differ- 
ent cooling curves for the various pairing types in Fig. 1221 
For singlet-state pairing, we see that neutrino emission is 
more strongly suppressed and the resulting temperature 
is slightly higher than the other types of pairing. This 
result confirms that the cooling behavior is not strongly 
dependent on the type of pairing of the baryons. 



17 



By comparing Fig. [20] and Fig. [TTl we see that the lu- 
minosity curves in the superfluid case are much closer 
to observations than in the non-superfluid case. How- 
ever, the curves do not agree with all the observational 
data except for the younger stars, namely number 1-5. 
This mismatch may lie in some important mechanisms 
which are neglected in the present calculation. For ex- 
ample, the magnetic field of the neutron star can affect 
the photon radiation at the surface [l^l such that it no 
longer obeys the standard black-body radiation law and 
shift the photon-radiation era to a later time [IHl- This 
reduction of photon radiation also decreases the luminos- 
ity, which is unfavorable according to the cooling curves 
shown in Fig. [20l However, there are some known heat- 
ing mechanisms which are not considered here, for exam- 
ple, due to magnetic field energy and rotational energy 
(for more details, see Ref. [l3] and references therein), 
such that the temperature and the luminosity can re- 
main higher. Besides, note that sources such as number 
10 and 11 may be old magnetars. As a consequence, the 
age estimates may not be correct, and the cooling his- 
tory may be anomalous. Both are a result of the decay 
of the strong magnetic field: age estimates assume mag- 
netic braking of rotation with a constant magnetic field, 
whereas the decay of the magnetic field results in heating 
of the neutron star [s^ . We discuss these issues in future 
work. 

Now we discuss the second ansatz, Tc = const. Since 
g-*(r) > jf ggi- ^ j.(o) = 10^ K we wiU have 
a lower Tc than in the previous ansatz and the baryons 
will not be superfluid everywhere inside the star ai t = Q. 
To make a proper comparison between the two, we set 
Tc = 1.9 X 10^ K, where the extra factor 1.9 is the middle 
value of e"**^'') inside the star (the outer value is taken at 
r = 7 km rather than at the surface since the direct Urea 
process starts around here, see Figs. l9l and ITO)) . With 
this value of Tc, at i = 0, the baryons near the sur- 
face are superfluid while those in the core are not, since 
is bigger in the core. Nevertheless we 
expect the integrated results will be comparable to those 
of the previous ansatz. As a comparison, in Fig. I23[ Cy 
and Qi, are shown for these two ansatze. We see that 
the difference between the two increases with decreas- 
ing temperature and only becomes significant when their 
values are very small. Therefore, it is not surprising to 
find that their cooling curves are similar, as is shown in 
Fig. [24l The two curves are quite close except that the 
second ansatz has a little lower luminosity, because with 
Tc = const the higher temperature in the core delays the 
suppression of the direct Urea process, such that the to- 
tal cooling rate is larger in the beginning. Furthermore, 
the result of the ansatz Tc = 10^ K mentioned above is 
also shown, where we can see that the high cooling rate 
due to the direct Urea process remains for a longer time 
and the luminosity is even more reduced after the neu- 
tron star has become completely superfiuid. Here we see 
again that different initial conditions have little influence 
on the later-stage behavior. 
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FIG. 23: The temperature dependence of Cv (solid curves for 
the singlet-state pairing contribution and dashed curves for 
the triplet-state pairing contribution) and Qv (dotted curves) 
calculated with the two different ansatze for Tc in the neutron 
star with the hadronic equation of state. For each quantity, 
the thick higher curve is for Tc = const and the thin lower 
curve is for Tce*^""^ = const, respectively. 
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FIG. 24: The luminosity calculated with two different ansatze 
for Tc in the neutron star with the hadronic equation of state. 
The solid curve is for Tce*'""^ = 10^ K and the dashed curve is 
for Tc = 1.9 X 10^ K, respectively. For comparison, the result 
with the ansatz Tc = 10® K is shown by the dotted curve. 
The observational data are again shown. 



Finally, we briefly present the results for the neutron 
stars with the mixed and nuclear equation of state. For 
the nuclear case, we always take the neutrons and pro- 
tons to pair in the triplet channel. The general argu- 
ments still apply, i.e., we only need to consider the direct 
Urea process and the photon radiation as the dominant 
cooling mechanisms. Without too much difference, the 
simple ansatz Tce*'-''^ = const will be used. The cool- 
ing behavior is mainly determined by the heat capacity, 
the direct Urea neutrino emissivity, and the photon ra- 
diation. These three parameters are not quite different 
among the various equations of state, even including the 
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FIG. 25: The luminosity of the neutron stars with the vari- 
ous equations of state with the effect of superfluidity included. 
The solid, dashed, and dotted curves represent the hadronic, 
mixed, and nuclear equation of state, respectively. The ob- 
servational data are again presented. 

effect of superfluidity. The dominant contribution to the 
heat capacity is due to the leptons. The direct Urea 
process is dominated by the triplet-state paired baryons 
which make up the majority of the star and are less af- 
fected by superfluidity. The photon radiation is again 
just black-body radiation. Therefore, we expect that the 
cooling process should be similar in these three cases, 
whose luminosity curves are shown in Fig. 1251 Compared 
to Fig.[T71 besides the overall increase in the luminosity as 
previously discussed, there is some difference in the order 
of magnitude of the luminosity among these three equa- 
tions of state. With superfluidity, the nuclear case always 
has the largest luminosity, which is due to the higher lep- 
ton fraction in the star. We find that the heat capacity of 
the leptons in the neutron star with the nuclear equation 
of state is about twice as large as that for the star with 
the hadronic or mixed equation of state. This larger lep- 
ton heat capacity becomes significant when the baryons 
(and the quarks) are highly superfluid. As can be seen in 
Eq. ([25)) . the larger the heat capacity is, the higher the 
corresponding temperature should be. This explains the 
largest luminosity of the nuclear case in Fig. 1251 Another 
difference is that, with superfluidity, the luminosity of the 
mixed case is always smaller than that of the hadronic 
case, while for the non-superfiuid case as in Fig.[T71 there 
is no such clear order. The reason also lies in the effect 
of superfluidity on the heat capacity, since in the mixed 
case the quark contribution is strongly suppressed due 
to superfluidity and completely neglected during our cal- 
culation. Of course, the neutrino emissivity from quarks 
is also reduced, but the reduction in the heat capacity 
is relatively stronger. This is because, without super- 
fluidity, quarks contribute less to the neutrino emissivity 
than baryons, but their contribution to the heat capacity 
is comparable. 

V. CONCLUSION AND DISCUSSION 

We have carefully compared the thermal evolution 
of different types of neutron stars, namely with the 



hadronic, the mixed phase of hadronic and strange quark 
matter, and the nuclear equation of state. We find that 
the direct Urea process is open in all of these cases and 
thus results in relatively fast cooling behavior. Although 
the details concerning the heat capacity and neutrino 
emissivity can be rather different in these cases, the cool- 
ing curves are quite similar after the stars become isother- 
mal. However, the behavior in the early stages before the 
stars become isothermal could be significantly different, 
but such a study requires the knowledge of the thermal 
conductivity of these complex systems. The geometrical 
structure of the mixed phase is also expected to play an 
important role, although no decisive conclusion has yet 
been drawn. 

As we have seen, the fast cooling in the non-superfiuid 
case did not agree with observations. In order to rem- 
edy this discrepancy, superfiuidity was introduced, which 
significantly reduces the efficiency of the direct Urea pro- 
cess as well as the heat capacity. The resulting cooling 
curve is much closer to the observational data. We also 
found that the particular pairing type of the superfluid 
baryons is not very important to the thermal evolution 
of the star. The thermal evolution after the star becomes 
isothermal is not strongly dependent on the initial tem- 
perature of the star or on the critical temperature related 
to superfluidity, as long as they are within reasonable 
ranges. The robustness of the results with superfiuidity 
is quite helpful in order to remove the uncertainties con- 
cerning baryon superfluidity at high density. Note that 
the cooling process is almost completely determined by 
the direct Urea process and photon radiation even after 
including the effects due to superfluidity. The cooling 
curves of the neutron stars with the three equations of 
state are still quite similar when superfluidity is included. 
We expect that a calculation including the magnetic field 
and rotation gives an even better agreement with the 
observational data. Besides, since our nuclear model is 
geared originally towards nuclear matter near the satu- 
ration point, a further improvement of it may play an 
important role in getting a better agreement with the 
observations. By incorporating the properties at higher 
densities, e.g., the coupling constants for hyperons, it is 
possible to get cooling curves covering most of the obser- 
vational data, as in Ref. [l^ . 
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